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Strongly interacting systems have been conjectured to spontaneously develop current carrying 
ground states under certain conditions. We conclusively demonstrate the existence of a commensu- 
rate staggered interlayer current phase in a bi-layer model by using the recently discovered quantum 
Monte-Carlo algorithm without the sign problem. A pseudospin SU(2) algebra and the corre- 
sponding anisotropic spin-1 Heisenberg model are constructed to show the competition among the 
staggered interlayer current, rung singlet and charge density wave phases. 

PACS numbers: 71.10.Fd,71.10.Hf, 71.30.+h, 74.20. Mn 



Strongly correlated systems can spontaneously break 
symmetries of the microscopic Hamiltonian. A partic- 
ularly interesting class of ground states spontaneously 
break the time reversal symmetry and carry a persistent 
current in the ground state. Such states are known by 
different synonyms, e.g. the orbital antiferromagnetic 
phase (OAF), the staggered flux (SF) or the D-density 
wave (DDW) phase. In the context of high Tc super- 
conductivity, these current carrying ground states have 
been proposed as competing states for the pseudogap 
phase[l S H II H @ The SF or the DDW phase has 
the attractive feature that the nodal quasi-particles have 
an energy spectrum similar to that of the d— wave super- 
conducting state. 

Whenever new ground states are proposed, it is impor- 
tant to establish for which microscopic Hamiltonian such 
states are realized. Because of their relative simplicity 
and availability of reliable analytical and numerical meth- 
ods, the ladder system has been used as a theoretical lab- 
oratory to investigate the DDW phase. Weak coupling 
bosonization methods combined with the renormaliza- 
tion group (RG) analysis on extended two-leg Hubbard 
ladders show the existence of commensurate DDW phase 
at half-filling 0, 0, and incommensurate power law 
fluctuating DDW order away from half-fillingp. Hoi . 
While the DDW state does not appear to be the ground 
state of the t-J ladder H^j ]. numerical works using 
the density matrix renormalization group(DMRG) found 
commensurate DDW order at half-filling l4 and incom- 
mensurate DDW order at low doping |l5| in a ladder 
model first proposed by Scalapino, Zhang and Hanke[l^|. 
The work of Schollwock et al has generated significant 
interest in connection with the DDW proposal for the 
cupratesQ. 

To the best of our knowledge, the existence of a current 
carrying ground state has not been conclusively demon- 
strated in any higher dimensional models. Following the 



FIG. 1: (a) Sketch of a staggered interlayer current (SIC) 
phase. For clarity, we do not show the bottom layer current. 
By conservation, each site acts as a source or drain for the cur- 
rent within the bi-layers. (b) Top view of the the bi-layer. (c) 
Sketch of the SF or the DDW current pattern for comparison. 



insights we learned from the ID systems, we investigate 
the current carrying ground state in a bi-layer version of 
the model constructed in Ref. . This model was orig- 
inally constructed and extensively investigated because 
of the exact 5*0(5) symmetry when coupling constants 
satisfy a simple relation, and is commonly referred to 
as the SZH model [Tl El E1E3- Here we show that 
the recently discovered fcrmionic quantum Monte Carlo 
(QMC) algorithm without the sign problem[2l| can also 
be applied to this model at and away from half-filling, for 
a large set of parameters, including purely repulsive in- 
teractions. Using this highly accurate numerical method, 
we can conclusively demonstrate the existence of a cur- 
rent carrying ground state in this model. The current 
carrying ground state is illustrated in Fig. ^ with stag- 
gered interlayer currents (SIC) between the bi-layers and 
alternating source to drain currents within the bi-layers. 
Viewed from the top of the bi-layers, this current pat- 
tern is different from the SF or the DDW current pat- 
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tern, since it has a s-wave symmetry. While the SF or 
the DDW currents are divergence free within the layer, 
the SIC current is curl free within the layer. These two 
flow patterns can be considered as dual to each other in 
two dimensions. In this paper, we shall first discuss the 



physics of the SIC phase by mapping onto an effective 
spin one Heisenberg model, and then proceed with the 
QMC results. 

The Hamiltonian for the SZH model 0] generalized 
straightforwardly to the bi-layer system reads 



H = -t\\ { c \a C j° + 4a d 0,o + h - c -} ~ *± { C i,<A* + h - C } ~ M \ C \<r Ci <° + (/ - } + J J2^' C ' 

(ij) i i i 

+ Uj2( n i,hc - l/2)K,|,c - 1/2) + (n»,T,d - l/2)(n u , d - 1/2) + F^(n i>c - l)Kd - 1), (1) 

i i 

i 



where c and d denotes fermionic operators in the upper 
and the lower layers, respectively, a corresponds to up 
and down spins. At half-filling, fi — 0, and the model 
is particle-hole symmetric, — 1 sets the unit of en- 
ergy. Before discussing the SIC phase, we first discuss 
some general properties of the SZH model which were 
not known before. The SZH model was known to have a 
50(5) symmetry when J = 4(Z7 + V) and /i = 0, which 
unifies antiferromagnetism with superconductivity ^(|. 
Remarkably, it also has another SO (5) symmetry when 

J = A(U-V), t ± = 0, (2) 

valid for all filling factors. In order to distinguish be- 
tween the two different 50(5) symmetries, we call the 
former the particle-particle 5*0(5) symmetry, denoted by 
SO(5) pp , and call the latter the particle-hole 5*0(5) sym- 
metry, denoted by SO(5) p h- The mathematical structure 
associated with the SO(5) p h algebra, not necessarily the 
symmetry itself, plays a crucial role in constructing the 
fermionic QMC algorithm without minus sign problem. 

We construct a four component fermion field \& = 
jcn-, d„} . Using the five Dirac T a matrices given in 
Ref.[21|, we construct the fermion bi-linears 

na = * t y* L ab =¥^ (3) 

It is straightforward to check that [H, L ab \ — when Eq. 
@ is satisfied, thus demonstrating the exact SO(5) p h 
symmetry. The SZH model can be mapped exactly to the 
spin 3/2 Hubbard modelj^, by the identification Cf = 
C3/2, C[ = C1/2, d| = c_ 1/2 , di = c_ 3 /2, and the SO(5) ph 
symmetry maps exactly onto the 50(5) symmetry of the 
spin 3/2 Hubbard model. Because of the exact mapping 
from the SZH model to the spin 3/2 Hubbard model, we 
are able to use the QMC algorithm discovered in Ref. 
|2l|. which works without the minus sign problem in a 
large parameter regime. 

In studying the strong coupling phase diagram, SZH 
identified the Eq phase where the rung singlet state, de- 
picted in Fig. is the lowest energy state, and the E3 
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FIG. 2: The double occupancy state a) and c) and the rung 
singlet (b). a), b), c) are spin SU(2) singlets and form the 
triplet representation of the pseudospin SU(2) group. 

phase, where the CDW states, depicted in Fig. 01 and 
are the lowest energy states. The new insight gained 
from Ref. 0, reveals that the competition between 
these two phases could result in the DDW phase. In view 
of this insight, let us consider the following operators 

ni {i) = i/2j2H(i)d a (i)-dl(i)c a (i)}, 

(J 

n 6 (<) = 1/2 J2H(i)d*(i)+ 4 
Q(i) = L 15 = l/2^{4(i)c CT (i)-4(iK(i)}, 

a 

where a are the Pauli matrices. These operators de- 
scribe rung current [ri\), rung kinetic energy (715) and 
the CDW order parameter iff). These three operators 
form a pseudo-spin SU(2) algebra which are important 
for our discussion of the SIC phase. 

There are 16 states on each rung, including 8 bosonic 
states with particle number 0, 2 or 4 and eight fermionic 
states with particle number 1 or 3. We are interested 
in the three rung states shown in Fig. EjL,b,c, which 
form a spin-1 representation of the pseudospin SU(2) 
algebra defined above, with eigenvalues Q = 1,0,-1. 
n\ ± in.5 act as pseudospin raising and lowering operators 
which connect these three states to each other. At half- 
filling and under the condition that max({7, V — 3/4 J) < 
min(y + J/4, U + 2V, U/2 + V), these are the three low- 
est energy states. Furthermore, at U = V — 3/4 J ', these 
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FIG. 3: Phase diagram in the strong coupling limit. Both 
S0(5) lines are shown as well as QMC region with no minus- 
sign problem for any filling (hatched area): g > 0, g > 
and U c > 0. There is also another region with V < (not 
shown). In the yellow region, the low-energy bosonic states 
are a, b and c shown in Fig. |21 This is where we expect 
the competition between the SIC and the rung-singlet phase. 
Black dots correspond to models for which we have performed 
QMC simulations. 

three states are degenerate. In the strong coupling limit, 
we can construct an effective theory to describe the low 
energy physics by using a pseudospin-1 antiferromagnetic 
Heiscnbcrg model 

Hex — J P ^2 { n 5(i) n 5(j) + ni(i)ni 

(j)+Q(i)Q(j)},(4) 

with J p — 2t 2 /(V + | J). Several terms break the pseu- 
dospin SU(2) symmetry. The intra-rung hopping ij_ term 
acts as an uniform external magnetic field which couples 
to 715. Also, the deviation of U from V — 3/4 J removes 
the degeneracy between a), c) and b) states. These sym- 
metry breaking terms are described by the on-site part 
as 

H on = ]T{-2^n 5 (i,) + A[/(Q 2 (z)-l/2)} (5) 

i 

where AU = U — (V — 3/4 J). The nonzero value of 
AU also gives different corrections to the three exchange 
terms at the order of J p AU/U. We will neglect these 
corrections below because the more important symmetry 
breaking effect from AU has already been taking into 
account in the on-site part. H = H ex + H on describes 
a 2D antiferromagnetic spin one Heisenberg model in an 
uniform magnetic field t±, with either easy axis (AU < 0) 
or easy plane (AU > 0) anisotropy. 

For the easy axis case with AU < 0, the effective 
Hamiltonian reduces to an Ising model with Q = ±1 
states, in a transverse magnetic field IsjUlif. For t± = 0, 
and AU > 0, the rung singlet state (b) has the lowest 



energy. However, in this case, there is a competition be- 
tween the AU > term and the Heisenberg exchange 
term J p . For AU > zJ p , where z = 4 is the coordination 
number, the ground state is a featureless Mott insulating 
state which can be described as a product of the rung sin- 
glet state on each site. On the other hand, for AU < zJ p , 
it is more favorable to form linear combinations between 
the (a), (b) and (c) states, such that a staggered ground 
state expectation value of (n\) and (n§) is spontaneously 
developed, thus lowering the Heisenberg exchange energy 
H ex . In this case, and for t± = 0, the pseudo-spin vector 
can lie along in any direction in the (711,715) plane. At 
AU = 0, a finite value of t± > corresponds to a pseudo- 
spin magnetic field along the n 5 direction, which creates 
an easy plane in the (ni,Q) space. The antiferromag- 
netic component of the pseudo-spin moment lies in the 
(ni , Q) plane, but the uniform component of the pseudo- 
spin moment points along the 715 direction. The pseudo- 
spin moment becomes fully polarized when t± > §Jp, 
and the antiferromagnetic component vanishes beyond 
this point. We see that t± > favors the (n\,Q) easy 
plane while AU < zJ p favors the (711,715) easy plane, 
therefore, when both conditions are satisfied, the inter- 
section between the two easy planes, namely the n\ easy 
axis, is selected. This is exactly the staggered inter-layer 
current (SIC) order. Combining all these considerations, 
we can summarize the subtle criteria for the SIC phase 
as 

V - < U < mm(V + j,2V), V > 

to. < \zJ p ^l -(AU/z J p ) 2 , AU < z,J p , (6) 

The first two robust conditions ensure that the (a), (b) 
and (c) states are the lowest and next lowest energy states 
among the 16 states on the rung, while the last two condi- 
tions are the rough mean field estimates discussed above. 

On Fig. |21 we show some specific regions on the phase 
diagram, obtained in the strong coupling limit. There 
are two additional axis for and tj_. If t|| and/or tj_ 
gets larger, we can expect some phases to have larger or 
smaller extension. In the case of ladders, a similar phase 
diagram has been proposed 0,^3- I n or der to obtain sig- 
nificant current correlations, one should be close enough 
from the line V = U + 3/4 J shown on Fig. |31 where states 
a, b and c become degenerate. 

Now we proceed to discuss the QMC calculation of the 
SIC phase. We first express the interaction terms of the 
SZH model as 

H in t = -g(n\ + nl)-g\nl+nl + nl)-U c (n-2) 2 , (7) 

up to a constant term. Here 4U C = — U — 3V + 3J/4, 
4g = V - U + 3 J/4 and 4g' = U - V + J/4. The SO(5) ph 
symmetry is clearly recovered when g — g' , i.e., when 
U = V + J/4. We now introduce auxiliary Hubbard- 
Stratonovich fields to decouple each of the three terms 
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o-o J(Q)/N 
b-b J(L/2,L/2) 
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FIG. 4: Parameters are t± = 0.1, f7 = 0, V = 0.5, J = 2.0 
and correspond to g = 0.25, </ = and U c = 0. Scaling of 
J7(Q)/N and J(L/2, L/2) vs 1/L showing almost no finite-size 
effects and proving long-range order in the thermodynamic 
limit. (Inset shows the convergence of J(Q)/N with the pro- 
jection parameter 6. Typically the GS value is obtained for 
9 = 20). 



above. Wu, Hu and Zhang|2lJ have shown that the QMC 
algorithm is free of the minus sign problem provided all 
three coefficients, g, g' and U c are positive, ft corre- 
sponds to a wedge in the phase diagram shown on Fig. |3 
and most remarkably, it includes a region with purely re- 
pulsive interactions, where U, V and J are all positive. 
A simpler case containing only n\ interaction, which ex- 
plicitly breaks the SU(2) spin rotation invariance, has 
been studied in another context j^. The ground-state 
(GS) properties of our model are conveniently studied 
with the projector auxiliary field QMC algorithm. The 
basic idea is to apply the operator exp(—9H) to a trial 
state. When 9 becomes large enough and with a proper 
normalization, this state converges exponentially to the 
GS. Details of the algorithm may be found in |23J]. The 
Trotter discretization was chosen to be At = 0.1 but we 
checked that it does not change the results. 

We compute correlations between rung currents n± (f) 
and perform its Fourier transform 



(8) 



The strongest signal in the Fourier transform is found 
for Q = (n,w), suggesting a staggered current pattern 
as shown in Fig. ^ This quantity converges to its GS 
value as the projector parameter 9 increases as shown in 
the inset of Fig. fn order to obtain information in the 
thermodynamic limit, one has to make an extrapolation 
of these GS values with a 1/1 finite-size scaling, where 
L is the linear size (L = 4, 6 and 8 in our simulations). 

Following our previous mean-field arguments, in order 
to prefer a phase with staggered current, we choose g > g' 
and U c = 0, with a small t±. As shown on Fig.0]for U = 
0, V = 0.5 and J — 2 when t± is small at 0.1, our values 
are rather constant with size, as expected in an Ising-like 
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FIG. 5: Finite-size scaling of the current correlations 
J{Q)/N showing no long-range order in the thermodynamic 
limit. Parameters are : (i) same as Fig.^Jexcept for t± = 0.5 
(ii) U = V = 0.3, J = 1.6 and t± = 0.5 at half-filling; (iii) 
U = 0.75, V = 0, J = 1 and t± = at 1/8-doping. Typically 
the GS value is obtained for 6 — 20. 



phase. Both the largest distance real-space correlations 
J(L/2, L/2) and the Fourier transform J(Q)/N converge 
to the same finite value (within our error bars) , meaning 
long-range order in the thermodynamic limit. 

As expected from our analytical estimates in ©, if 
AU or t± gets too large, long-range order disappear as 
shown on Fig. [SJ Since we can also perform the QMC 
simulation at finite-doping without the sign problem, we 
have chosen to work at 1/8-doping for some parameters 
shown on Fig.|SJ Again, rung-current correlations vanish 
in the thermodynamic limit since the Fermi surface is not 
nested anymore. From the analytical estimates based on 
the mapping to the spin one antiferromagnetic Heisen- 
berg model and the detailed QMC calculations shown 
above, we can conclusively demonstrate the existence of 
the SIC phase at half-filling, and also note that this is a 
rather subtle phase which can be easily destabilized by 
large U and doping. 
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